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Summary 

The theory and user instructions for an optimization code 
based on the method of feasible directions are presented. The 
code was written for wide distribution and ease of attachment 
to other simulation software. Although the theory of the method 
of feasible direction was developed in the 1 960’ s, many consid- 
erations are involved in its actual implementation as a computer 
code. Included in the code are a number of features to improve 
robustness in optimization. The search direction is obtained by 
solving a quadratic program using an interior method based on 
Karmarkar’s algorithm. The theory is discussed, focusing on 
the important and often overlooked role played by the various 
parameters guiding the iterations within the program. Also 
discussed is a robust approach for handling infeasible starting 
points. Thecode was validated by solving a variety of structural 
optimization test problems that have known solutions obtained 
by other optimization codes. It has been observed that this code 
is accurate and robust: it has solved a variety of problems from 
different starting points. However, the code is inefficient in that 
it takes considerable CPU time as compared with certain other 
available codes. Further work is required to improve its effi- 
ciency while retaining its robustness. 

Introduction 

The theoretical basis for the method of feasible directions 
(MFD) was originally developed by Zoutendyk (ref. 1). In 
engineering, the text by Vanderplaats (ref. 2) contains details 
on the actual implementation of the method. Vanderplaats’ 
algorithm, based on the method of feasible directions and 
published in 1983 (ref. 3) has had a great impact on engineering 
optimization. The program, as discussed here, is based on the 
same basic theory as that cited in references 2 
and 3 although the implementation aspects are different. A 


major challenge in developing software for engineering opti- 
mization is in implementing conceptual algorithms. The main 
purpose of this report is to describe in detail various implemen- 
tation aspects. It also serves as a users guide and describes the 
role played by the various parameters guiding the iterations in 
the program. 

The basic steps in the method of feasible directions involve 
solving a quadratic program (QP) to find the direction vector 
and then to find the step size along this direction by performing 
a constrained one-dimensional line search. In the implementa- 
tion herein, an interior method based on Karmarkar’s approach 
(ref.4) is used to solve the QP problem for direction finding. 
Interior methods have been found to be effective for problems 
with a large number of variables. Also, special attention is 
given to the treatment of infeasible starting designs, to the 
dynamic adjustment of the constraint thickness parameter 
during the iterations, and to line search strategies. The optimi- 
zation problem follows: 


Minimize 

F(DV) 

(1) 

subject to 

Gy(DV) < 0 7 = 1,..., NCON 

(2) 

and 

DVLj < DVj < DVUj i = 1,..., NDV 

(3) 


where F is an objective function to be minimized subject to the 
constraints Gj < 0, DV is an (NDVX1) vector of design 
variables, NDV is the number of design variables, NCON is the 
number of constraints, and DVL and DVU are the lower and 
upper bounds of the design variables, respectively. Equality 
constraints may be included by means of a penalty function 
(ref. 2). 

Here, the capital variable names DV, DVL, and DVU are 
those that are actually used in the program. The default values 
for the various parameters are given in appendix A. 
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Direction Finding and Line Search From 
Feasible Design Points 

Direction Vector 

Let DV° be a current design point which satisfies con- 
straints (2) and (3). The direction finding problem is to find a 
search direction d that will ( 1 ) point into the feasible region, or 
is “a feasible direction,” and (2) reduce the objective function, 
or is “usable.” If VF • d < 0 and VGj * d < 0 for each active 
constraint, then the direction vector d is usable-feasible. To find 
a direction that is both feasible and usable, the following 
subproblem is posed: 


Minimize 

e 



subject to 

VF 

d < 6 



VGj 

• d <Pj 0, j in J [./, 

(4) 


1/2 

VI 

E-n 

T3 



In problem (4 \Jep is an active set defined by the union of active 
constraints and active bounds 


where A is a matrix of dimension (NAC+1 X NDV+1), NAC 
equals the total number of active constraints (including active 
bounds). All gradient vectors in (6) have been normalized to be 
unit vectors. The problem in (4) can now be written as follows: 


Minimize 

p r y 


subject to 

A r y <0 

(7) 

and 

1/2 y T y < 1 



The dual formulation corresponding to problem (7) can be 
written as (ref 2) 

Minimize M2 pJ [A^Aj p + p^A^pj 

( 8 ) 

subject to p > 0 

where p is an (NAC+1 X 1) vector. Problem (8) is a QP 
subproblem which can be solved in a number of ways. Once the 
QP is solved, the solution to (6) is obtained from 

y = -p - Ap (9) 

which gives the search direction d. 

Solution of the QP Using an Interior Approach 


J EP - {j : Gj + EP > oj u 
{Active lower/upper bounds} (5) 

where 0 is an artificial variable, pj is the push-off factor 
associated with the y th constraint, and EP is the constraint 
thickness parameter. Again, see appendix A for the default 
value of EP and the values of the other parameters. The 
following notation is now introduced: 


y = [d, 0] T = an (NDV + 1x1) vector 

p = [0, ..., 0 , \} T = an (NDV + 1x1) vector 

‘ VC,, -ft 

VG 2 , ~~p2 


a t = 


V^NAO -&NAC 


VF, -1 


( 6 ) 


Interior methods gained greater acceptance after Karmarkar’ s 
paper in 1984. A variation of Karmarkar’ s approach, called the 
“linear affine scaling” algorithm, has also been developed and 
applied to linear programming problems (ref. 4). The approach 
was extended here to solve the QP in (8). The algorithm begins 
by first setting an initial estimate, p° > 0 . The main steps are 

( 1 ) A transformation to a new set of variables p new is defined 
by 

p = D p new 

where D is a diagonal matrix whose elements are equal to the 
components of p° = current value of p. Thus, the current value 
of the new variables will always equal [1,...,1] . 

(2) The search direction s is defined to be the steepest descent 
direction in the new space of variables 

s = -D (A r Ap + A r pj 

(3) The step size p along s is obtained from 

p = Minimum (pj, p 2 ) 
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where p\ is 98 percent of the maximum step to the boundary (to 
ensure that ft is positive), and p 2 is the minimum of the 
quadratic function along s. The variables at k th iteration are 
updated from 

\X k+I = fi k + p D s (10) 

and we go back to step 1 . These steps are repeated until changes 
in the objective function of the QP are relatively small for three 
consecutive iterations. 

Step Size Determination 

At a feasible point DV°, we have now found a search 
direction d which is both usable and feasible. The step size a 
along d is determined in the program as follows. First, an initial 
step ai n j t > a mjn > 0 is chosen based on the following criteria: 

a init = Minimum (a B , a F , a ML ) (11) 

where 

a ff maximum step based on lower and upper limits on DV 

ap value of a that reduces the objective by 5 percent 

a ML value of a that results in a 1 5-percent maximum change 

in a design variable 

Once the initial step is determined, a move is initiated based on 
the golden section ratio GR of 0.61 8. Thus, steps ai n j t , a m \ { /GR, 
ainit/(G/?) 2 , ... are taken along d, with a < ag. At every step in 
this move, the design variables are updated as 

DV(a) = DV° + a • d (12) 

and the objective and constraints are evaluated. If at any stage 
the constraints are violated, then the previous feasible point and 
the current infeasible point define an interval containing the 
zero of the constraint function. A bisection scheme is used to 
determine this zero, which determines the maximum limit on 
the step size. The bisection process is terminated when the 
maximum constraint value G max satisfies 

— ^low — ^max — ^high 03) 

where Ghigh is set equal to zero and should not be changed ; Gj ow 
is set equal to 0.001 by default. If the objective function is 
linear, then the root of the constraint function is the desired step 
size «o- 

The other possibility is that the objective function is nonlin- 
ear along d. Thus, we also monitor the directional derivative of 


the objective function along d at every step during the move 
from 

F{a) = VF[DV(a)]d 

Since d is a usable direction, F < 0 at a = 0. If F > 0 at any 
feasible a during the move, then the previous point where F is 
negative and the current point where F is positive defines an 
interval containing the minimum of F. A bisection scheme is 
used based on computing F{a) to determine the minimum. The 
bisection scheme is terminated when the interval length is less 
than a small number (TOLINT). 

Once the step size a 0 is determined, the new design variables 
are obtained from 

DV 1 =DV°+a 0 d 

We reset DV 1 to equal DV° and repeat the process. That is, we 
again find a new search direction and step size. 

Selection and Adjustment of Push-Off Factors 

The push-off factors fy entering into the direction finding 
problem in (4) are determined from equation (2): 

Pj=p 0 (l + Gj/EP ) 2 (14) 

where j3q is equal to 1 by default. The push-off factors are set 
equal to zero for bounds and linear constraints. Push-off factors 
are also adjusted based on the following. 

The solution of the QP in (8) and in equation (9) gives us a 
vector d. However, if VF d < 0 but VG p d> 0 for some active 
p, then the search direction is usable but not feasible. In this 
case, the corresponding push-off factor fi p is increased a little 
and d recomputed. This increase is only attempted five times 
after which EP is reduced as discussed in the following section. 
In theory, the vector d from the QP solution should equal zero 
in such situations where a usable-feasible vector cannot be 
determined. However, it should be noted that thfc interior 
method yields solutions which are interior to the boundary. 
Round-off errors may also lead to such situations. 

Adjustment of Constraint Thickness Parameter EP 

The constraint thickness parameter EP is used to determine 
the active set as given in (5). The default value at the beginning 
of each design iteration is EP = EP 0 = 0.01 . If at any stage, the 
QP solution yields a search direction d that is not usable- 
feasible, then EP is reduced. This happens particularly as the 
optimum is approached. We iteratively reduce EP as 

EP new = EP old /3 (15) 


3 



which is done at most five times after which the program 
terminates with the final solution. 


expressed in the form of a QP as given in (8) with only one 
difference: the last row in the A r matrix in (6) is deleted. 


Stopping Criteria 


Line Search 


The three stopping criteria used are as follows: 

Kuhrt-Tucker condition. — First, determine if Idl < 
TOLKKT, where TOLKKT is a small number (=10"* by de- 
fault) and d is the search direction. Usually, this stop criterion 
is satisfied only for simple and small problems. 

Changes in the objective function . — Once we have ob- 
tained a usable-feasible direction, a line search is performed 
to obtain a new point with a lower value for the objective 
function. If either absolute changes in For relative changes 
in F are small for three consecutive iterations, then the pro- 
gram terminates. Thus, we check 


abs (change in F) < TOLABS 


abs (change in F/F) TOLREL 


(16) 


Assume that the QP in (17) has been solved and a search 
direction d has been obtained for constraint correction. First, 
is determined along d so that lower and upper bounds on the 
design variables are satisfied. Within the interval [0, a#], a 
move is initiated to determine the step size a$. At any point in 
this move, if the maximum constraint value G max < 0, then the 
constraint correction phase terminates and we switch over to 
the algorithm in the section Direction Vector. 

If the violation continues to reduce along d and a = a#, then 
we set a equal to a# and compute a new search direction for 
constraint correction. On the other hand, if we notice the 
violations to reduce and then increase, then the interval contain- 
ing ocq can be bracketed. That is, if there are three points a j, ai, 
c *3 for which G maxi > G maX2 > G maX3 , then we have a\ < a$ < 
a$. A golden section search is then performed within this 
interval to obtain the minimum of G max . A new search direction 
is again computed at the new point and this process is repeated 
until a feasible design is obtained. 


Default values of TOLABS and TOLREL are 1(T 6 . Of 
course, if F has a very large magnitude, then TOLREL can be 
set equal to a smaller value. 

Successive reductions in the constraint thickness param- 
eter : — If, after successive reductions in EP, there is no usable- 
feasible direction, then the program terminates with the final 
solution. 

Direction Finding and Line Search From 
an Infeasible Starting Point 

If the starting design DV° is infeasible, with at least one 
constraint Gy > 0, the violations are corrected first. Then 
proceed with the algorithm as described in the preceding 
section. At an infeasible point, the dot products of the violated 
constraint functions with the direction vector should be nega- 
tive. Thus, we define the following subproblem for direction 
finding: 

Minimize 0 


Push-Off Factors and Constraint Thickness Parameter 

For the constraint correction subproblem in equation (17), 
push-off factors fy are computed from equation ( 1 4). However, 
if any )3y > POM AX, then /3y is set equal to POMAX. The default 
value of POMAX is 50. 

Sometimes, it is not possible to find a direction vector that 
makes negative dot products with all the active or violated 
constraints. In this case, the constraint thickness parameter EP 
is decreased in order to reduce the number of constraints in the 
active set. For lower and upper bounds, EP is reduced as given 
in equation (15). For the constraints Gy, EP is reduced by the 
formula 

EP new =EP_0.5 (EP + C^) (18) 

Since G max > 0, EP can take on a negative value. Equation ( 1 8) 
is applied a maximum of five times after which an error 
message is printed out to the effect that no feasible design can 
be found. The strategy in equation (18) is new and is quite 
effective. 


subject to V Gy d < /3y 6 , j in J EP (1 7) 

1 / 2 d r d < 1 

In reference 2, the objective function also enters into the 
subproblem. Herein, the goal is to exclusively correct con- 
straint violations. As it turns out, the problem in (17) can be 


Users Guide 

The program MFD has been written in QB ASIC as well as in 
Fortran. This section describes the user-supplied subroutine 
with an example (source listing in appendixes B and C). 
Appendix A contains a description of the main parameters used 
by the program together with their default values which may be 
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changed in the user routine. 

User-Supplied Subroutine 

To run the program, three steps are involved: 

FIRST . — Upon initiation, the program will ask for the val- 
ues of NDV and NCON, which are the number of design 
variables and the number of constraints. 

SECOND — In SUBROUTINE USER, there are four “CALL 
XXXX” statements which call the user-supplied subroutine. 
The user should give the appropriate name of the user-supplied 
subroutine in these CALL statements. 

THIRD , — The user must supply a subroutine. This is now 
explained. 

The user routine will be called with INF equal to -1 , 1 , or 2. 

INF=-1 . — This is the first call to the user routine. The initial 
values of the variables together with their lower/upper bounds 
have to be defined in arrays DV, DVL, DVU, respectively. 
Also, the user may change the defaults of the parameters as 
discussed subsequently in appendix A. For example, IGRAD = 1 
may be set indicating that gradients will not be supplied, the 
number of iterations ITLIM may be set, or LINEARF = l may 
be set if the objective function is linear. The general appearance 
of statements will be 

IGRAD j 

ITLIM > = See appendix B for explanations and 
LINEARF J default values. 

DV = Initial guess of the design variables 

DV = Lower bounds 

DVV = Upper bounds 

INF = L — The objective and constraint functions have to be 
defined for the given values of the variables contained in the 
array DV. That is, we need to supply 

F 

G(l) 

G(2) 

G(NCON) = 

INF = 2 . — If IGRAD = 0 (default), then gradients of the 
objective function and of the active constraints have to be 
supplied in the matrices DF and AA, respectively. The active 
constraint numbers are supplied to the routine in the array 
IACT(1),...,IACT(NAC). If NAC = 0 at any iteration, then 
there are no active constraints and only the objective function 


gradient DF needs to be computed. 

For example, let NDV = 3, NCON = 10, NAC = 2, IACT (1) 
= 3, IACT(2) = 9. Then, we must supply 

DF(1) = , DF(2) = , DF(3) = 

AA( 1 , 1 ) = , AA(2, 1 ) = , AA(3, 1 ) = ,- Gradient of third 
constraint in first column of [AA] 

AA( 1,2) = , AA(2,2) = , AA(3 ,2) = ,- Gradient of ninth 
constraint in second column 

Example Problem 

Consider the Rosen-Suzuki problem: 

Minimize F = XI 2 + X2 2 +2 X3 2 -X4 2 -5 XI 
-5X2-21 X3 + 7X4 + 100 

subject to 

XI 2 +X2 2 +X3 2 +X4 2 +X1-X2 + X3 -X4-8<0 
XI 2 +2 X2 2 +X3 2 +2 X4 2 -XI — X4 — 10 < 0 

2 XI 2 + X2 + X3 2 + 2 XI - X2 - X4 - 5 < 0 

The solution to this is (0. 1 5, 1.0, 1.87,-1.14) with F opl equal 
to53.6.Userroutines for this problem withIGRAD= 1 and also 
a spring design problem with IGRAD = 0 are provided in 
appendix B. Note that the common blocks are required in the 
routine along with the double precision statement. 

Test Problems 

The program MFD has successfully solved a variety of test 
problems. The results for only three truss optimization prob- 
lems will be given here. The results tally with those obtained 
from using other optimization codes from a test bed that was 
developed at NASA Lewis Research Center (ref. 5). 

Ten-Bar Truss 

This problem, given in reference 6, is to minimize the weight 
of the 10-bar truss shown in figure 1 . The data are as follows: 
Young’s modulus E= 10 7 psi;p = 0.1 lb/in. 3 ; allowable stress 
|| = 25 000 psi for elements 1 to 8 and 10; <7 a p = 75 000 psi 
for element 9; NDV = 10; NCON =10; the lower bounds on 
each area = 0.1 in. 2 ; the upper limit = 20.0 in. 2 . In figure l,load 
P = 100 000 lb. The initial design variables were 4.0 in. 2 
(infeasible) with the initial weight F° = 1678.6 lb. 
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The optimum design is as follows: 

F ° pt = 1498.3 lb 

DV opt ={7.92, 0.1, 8.10, 3.90, 0.1, 0.1, 5.80, 
5.52, 3.67, 0.14} in. 2 

(The cross-sectional areas given above are rounded to two 
decimal places.) 

Active set: stresses in elements 1,2,3,4,6,7,8,10 

The results from other codes are given in table I. Program MFD 
solved this problem in 127 iterations with 1007 function evalu- 
ations. Because this number may seem large for a 10-bar truss, 
detailed discussions on this aspect are presented in the next 
section. 


Tapered 10-Bar Truss 

This problem is a variation of the previous problem (see 
fig. 2). There are constraints on the stress in each of the 10 
elements and also on the displacements in the Y-direction of 
nodes 3 and 4. The constraints are to be satisfied under two 
loading conditions (table II). Thus, NDV = 10 and NCON = 24; 
E = 10 7 ; p = 0.1 lb/in. 3 ; cr a jj = 10 000 psi in each of the 10 
elements; the lower bounds = 0.1 in. 2 ; the upper bounds = 
100 in. 2 . The initial design is (l.,...,l.) in 2 with F° = 146 lb 
(infeasible). The optimum solution is 

F opt = 3269.0 lb 

DV opt ={57.34, 18.66, 2.25, 7.41, 36.13, 3.44, 

28,79, 19.69, 16.20, 8.94} in 2 

Active set for load case 1: stresses in elements 1,7,8 and 
displacement at node 3 

Active set for load case 2: stress in element 6 

MFD took 121 iterations and 1091 function evaluations. Table I 
compares the performance of other codes for this problem. 

Sixty-Bar Trussed Ring 

The trussed ring in figure 3 consists of 60 elements and is 
subjected to 3 loading conditions (table III). The data are 
E= 10^ psi; p = 0.1 lb/in. 3 ; a a p = 10 000 psi for each element; 
node 4 has a displacement limit of 1.5 in. in the Y-direction; the 
lower and upper limits on areas are 0. 1 and 100.0 in. 2 , respec- 
tively. The 60 element areas are linked to the 25 design 
variables given in table IV. Here, NDV = 25 and NCON = 1 83. 


TABLE I — OPTIMUM WEIGHTS FROM VARIOUS CODES 
[Number of function evaluations shown in paretheses.] 


Method 

Truss optimization problem 


Ten-bar 

Tapered ten-bar 

Ring 


Optimum weight, F t lb 


Method of feasible direction (MFD) 

1500 

(1007) 

3269 

(1091) 

345 

(1836) 

Sequential unconstrained 

1503 

3270 

345 

minimizational techinque (SUMT) a 

(243) 

(667) 

(369) 

Feasible direction (FD) a 

1498 

(73) 

3424 b 

(77) 

343 

(78) 

Sequential quadratic programming (SQP) a 

(c) 

3271 

(187) 

344 

(782) 


Reference 5. 

b A different starting point was necessary. 

C A feasible design using the default parameters in the code could not be obtained. 
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tively. The 60 element areas are linked to the 25 design 
variables given in table IV. Here, NDV = 25 and NCON = 1 83. 
The initial design is (l.,...,l.) in.^ with an associated weight of 
250 lb. The optimum design is 

F opt = 345.0 lb 


{1.19 

2.17 

0.14 

2.18 

2.13 

0.68 

2.20 

2.20 

1.08 

2.26 

DV°P‘ = 2.18 

0.26 

2.33 

1.25 

1.22 

0.88 

0.88 

1.23 

1.12 

1.15 

1.27 

1.15 

0.88 

1.23 

1.26} 


Active set for load case 1 : stresses in elements 25, 49, and 
displacement at node 4 

Active set for load case 2: stresses in elements 5, 14, 20 

Active set for load case 3: stress in element 1 1 

The optimum was obtained in 247 iterations with 1 836 function 
evaluations. Again, see table I for the performance of other 
codes. 


Discussion 


TABLE II.— LOADING 
SPECIFICATIONS FOR 
TAPERED 10-BAR TRUSS 


Load 

condition 

Node 

Load direction 

X 

Y 

Load, P . lb 

Case ! 

2 

60 000 

120 000 


3 

60 000 

60 000 


4 

17 500 

12 500 


5 

17 500 

25 000 

Case 2 

2 

0 

-50 000 


3 



-25 000 


4 



-37 500 


5 



-75 000 
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Figure 3. — Sixty-bar trussed ring. 


In this report, the detailed implementation of a feasible 
directions code was explained with particular emphasis on line 
search, infeasible starting points, and the role of various param- 
eters used in the algorithm. 

The robustness of the program was very satisfactory for a 
variety of test problems, both with regard to correcting an 
infeasible starting point and to making subsequent progress 
towards an optimum solution. Its robustness may be attributed 
to the use of the interior method for direction finding and to the 
use of a bisection scheme in the line search. 

The following factors contribute to the inefficiency of the 
code: 

(1) A plot of truss volume (cost) versus the number of 
iterations for the 10-bar truss problem is shown in figure 4. We 
see that substantial reductions take place in the first 40 itera- 
tions while the remaining 120 iterations are needed only to 
achieve the convergence criterion of 1CT 6 . This asymptotical 
characteristic is typical of most mathematical programming 
codes. In an interactive setting, the user can terminate the 
optimization run as soon as the cost-versus-iteration curve 
flattens out. 

(2) In the final design, there are no constraint violations 
whatsoever because of the line search control described by 
equation (13), which is in contrast to certain other programs that 
yield final designs with small violations. 

(3) Push-off factors are user-defined. Their default values 
are unity. Their selections can be optimized for categories of 
problems, truss being one. Such studies have not been carried 
out here. 

(4) Efficiency is sacrificed when pure bisection is used. In the 


7 



Volume of truss, in/ 


TABLE III.— LOADING 
SPECIFICATIONS FOR 
60-BAR TRUSSED RING 


Load 

condition 

Node 

Load direction 

X 

Y 

Load, P, lb 

Case 1 

1 

-10000 

0 


7 

9000 

0 

Case 2 

15 

-8 000 

3 000 


18 

-8 000 

3 000 

Case 3 

22 

-20 000 

10 000 


TABLE IV. — GROUPING OF 
ELEMENTS FOR 60-BAR 
TRUSSED RING 


Design 

group 

Element 


Design 

group 

Element 

1 

49-60 


13 

12,24 

2 

1, 13 


14 

25, 37 

3 

2, 14 


15 

26, 38 

4 

3, 15 




5 

4, 16 


16 

27, 39 




17 

28,40 

6 

5, 17 


18 

29,41 

7 

6, 18 


19 

30, 42 

8 

7, 19 


20 

31,43 

9 

8, 20 




10 

9,21 


21 

32 ,44 




22 

33,45 

11 

10, 22 


23 

34,46 

12 

11,23 


24 

35,47 




25 

36, 48 



test problems, a line search within the feasible region has taken 
between 5 and 10 function evaluations (analyses) as determined 
by taking the ratio of the number of function evaluations to the 
number of iterations. There is latitude here to improve the 
efficiency of the algorithm by using hybrid approaches which 
combine sectioning with polynomial interpolations for line 
search. 

Conclusions 

The method of feasible directions incorporates Karmarkar’s 
algorithm for the generation of search directions and the bisec- 
tion technique for improvements in one-dimensional search. 

The algorithm provides a feasible optimum design for prob- 
lems with complex constraint space even when the initial design 
is infeasible. 

For selected problems the computer code may require more 
computations than those that may be needed by other nonlinear 
programming software. However, there is latitude to improve 
the computational efficiency of the algorithm. 

The computer code successfully solved several structural 
optimization problems. 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, November 2, 1993 
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Appendix A 

Program Parameters and Their Default Values 

This section describes various parameters used by the program and their default values. The user may change any of these values 
in the user-supplied subroutine when INF = -1 as explained in the section User-Supplied Subroutine. (The important parameters 
are indicated with an asterisk.) 

ALPMIN / 0.001 / — minimum step length used in line search procedures 

ALPQP / 0.98 / — a parameter used to come near the boundary in the interior approach to solve the QP 
AMLIMIT, OBJFAC, ALPMIN /0.15, 0.05, 0.001/ — used during line search 

jjj 

EP /0.01/ — constraint thickness parameter (base value) 

EPBD /0.01/ — constraint thickness parameter for lower and upper bounds 

*EPSPCT, EPSMIN /0.01 , 0.001/ — parameters relevant only when IGRAD = 1 for divided-difference gradients (The gradient f is 
evaluated by (f(x+e) -f(x)}/e, where e = maximum (EPSPCT X abs(x), EPSMIN).) 

GLOW, GHIGH /0.001 , 0/ — used during line search to determine the zero of a constraint function (GHIGH should not be changed 
from its zero value.) 

* IGRAD 101 — IGRAD =0 means that the user will provide gradients of objective and active constraint functions in the user-supplied 

subroutine (see next section). IGRAD = 1 means that gradients will not be provided by the user; they will be internally calculated 
in the program using divided differences. See EPSPCT and EPSMIN 

ITEPMAX / 5 / — maximum number of reductions in the constraint thickness parameter EP to obtain a usable-feasible direction or 
to obtain a feasible direction during constraint correction 

ITLIM /1 50/ — maximum number of iterations 

ITLIMINF / 30 / — maximum number of iterations to obtain a feasible design from an infeasible starting point 
ITPUSHMAX / 5 / — maximum number of iterations to increase the push-off factors to obtain a usable-feasible search direction 
*LINEARF / 0 / — LINEARF = 0 if objective function is nonlinear and = 1 if objective function is linear 
LINEAR(J) /0,...,0/ — = 0 if constraint function is nonlinear and = 1 if constraint is linear 

* MAX ACT / 50 / — maximum number of active constraints; used for dimensioning purposes 

POMAX / 50. / — maximum push-off factor used during constraint correction phase for infeasible starting designs 
TOLINT / 0.001 / — used as stopping criterion during line search (associated with nonlinear objective functions) 

TOLKKT / 0.0001 / — stopping criterion used for checking Kuhn-Tucker conditions 

*TOLABS, TOLREL /l.e-6, l.e-6/ — absolute and relative tolerances used for stopping criteria (see (16)) 

TOLABS 1 , TOLREL1 , ITMAX / 1 .E-6, 1 .E-6, 100 / — stopping parameters associated with the interior method for solving the QP 
for the direction vector 
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Appendix B 

Fortran User Subroutines with IGRAD = 1 and IGRAD = 0 

ROSEN-SUZUKI PROBLEM : IGRAD = 1 (divided-difference gradients) 

SUBROUTINE ROSEN 
IMPLICIT REAL* 8( A-H,0-Z) 

COMMON /MFD1/ INF, NDV, NCON, NFV, IGRAD, EP, NAC, EPSPCT, EPSMIN, 

$ MAXACT, EPBD 

COMMON /MFD2/ DV(100), IACT(50), LINEAR(IOOO), F, G(1000), 

$ DF(100), AA(101,51) 

COMMON /MFD6/ DVO(IOO), FI, FM, Gl(1000), D(100), DVL(IOO), 

$ DVU( 1 00), BETA( 1 000) 

IF (INF .EQ. -1) THEN 
C NDV = 4, NCON = 3 


IGRAD = 1 for divided-difference gradients 

DV(1) = 1 
DV(2) = 2 
DV(3) = 3 

DV(4) = 4 DV contains initial values of design variables 

D VL( 1 ) = - 1 0 
DVL(2) = -10 
DVL(3) = -10 

DVL(4) = -10 lower bounds 

DVU(l) = 10 
DVU(2) = 10 
DVU(3) = 10 

DVU(4) = 10 upper bounds 


RETURN 
END IF 
B1 =s DV(1) 

B2 = DV(2) 

B3 = DV(3) 

B4 = DV(4) 

NFV = NFV + 1 

Objective and constraint functions 


F = B1 ** 2 + B2 ** 2 + 2 * B3 ** 2 - B4 ** 2 - 5 * B1 - 5 * 

$ B2-21 * B3 + 7.*B4 + 100. 

G(l) = B1 ** 2 + B2 ** 2 + B3 ** 2 + B4 ** 2 + B1 - B2 + B3 - B4-8 
G(2) = B1 ** 2 + 2 * B2 ** 2 + B3 ** 2 + 2 * B4 ** 2 - Bl- B4 - 10 
G(3) = 2 * Bl ** 2 + B2 ** 2 + B3 ** 2 + 2 * Bl - B2 - B4 - 5 
RETURN 
END 
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Appendix C 

A Spring Design Problem : IGRAD = 0 (User-Supplied Gradients) 

SUBROUTINE SPRING 
IMPLICIT REAL*8(A-H,0-Z) 

COMMON /MFD1/ INF, NDV, NCON, NFV, IGRAD, EP, NAC, EPSPCT, EPSMIN, 

$ MAXACT, EPBD 

COMMON /MFD2/ DV(100), IACT(50), LINE AR( 1 000), F, G(1000), 

$ DF(100), A A( 10 1,51) 

COMMON /MFD6/ DV0( 1 00), F 1 , FM, G 1 ( 1 000), D( 1 00), DVL( 1 00), 

$ DVU( 1 00), BETA( 1 000) 

IF (INF .EQ. -1) THEN 

C SPRING PROBLEM — NDV =3: NCON = 4 

IGRAD = 0 by default 

LINEAR(4) = 1 Constraint #4 is linear 

DV(1) = 1 
DV(2) = 2 
DV(3) = 3 
D VL( 1 ) = .05 
DVL(2) = . 1 
DVL(3) = 1 
DVU(l) = 1 
DVU(2) = 5 
DVU(3) = 50 
RETURN 
END IF 
B1 =DV(1) 

B2 = DV(2) 

B3 = DV(3) 

IF (INF .EQ. 2) GOTO 701 

INF = 1 ( supply F and Gj) 

NFV = NFV + 1 

F = (B3 + 2) * B1 ** 2 * B2 

G(l) = 1 -B2**3*B3/71875 /B1 ** 4 

G(2) = (4*B2**2-B1 * B2) / 12566 / (B2 * B1 **3-Bl **4) 

G(2) = G(2)+ 1 / 5108/BI **2-1 
G(3)= 1 - 140.45 /B2** 2 / B3 
G(4) = (B 1 +B2)/ 1.5- 1 
RETURN 

Objective gradient in DF 

701 DF(1) = (B3 + 2) * 2 * B1 * B2 

DF(2) = (B3 + 2) * B1 ** 2 
DF(3) = B2 * B1 ** 2 

IF (NAC .EQ. 0)RETURN ... NAC = number of active constraints 
DO 10011= 1, NAC 

Constraint gradients placed in columns of AA- matrix using active set IACT 

IF (IACT(II) .EQ. 1 ) THEN 
AA(1, II) = 4 * B2 ** 3 * B3 / 71875 / B 1 ** 5 
AA(2, II) = -3 * B2 ** 2 * B3 / 71875 / B1 ** 4 
AA(3, II) = -B2 ** 3 / 7 1 875 / B 1 ** 4 
END IF 


n 



IF (IACT(n) .EQ. 2) THEN 
Cl = 12566 * (B2 * B1 ** 3 - B1 ** 4) 

C2 = 4 * B2 ** 2 - Bl * B2 
C3 = 5108 * B1 **2 

CC=-C1 *B2-C2* 12566 * (3 * B2 * B1 **2-4*Bl ** 3) 

CC = CC/C1 ** 2 

CC = CC - 2 / 5 108 / B 1 ** 3 

AA(1, II) = CC 

CC = (C1 * (8 * B2 - Bl) - C2 * 12566 *B1 **3)/Cl ** 2 
AA(2, n) = CC 
AA(3, II) = 0 
END IF 

IF (IACT(H) .EQ. 3) THEN 
AA(1, II) = -140.45 / B2 ** 2 / B3 
AA(2, II) = 2 * 140.45 *B1/B2** 3 /B3 
AA(3, II) = 140.45 * Bl / B2 ** 2 / B3 ** 2 
END IF 

IF (IACT(II) EQ. 4) THEN 
AA( 1 , II) = 1 / 1.5 
AA(2, II) = 1 / 1 .5 
AA(3, II) = 0 
END IF 
100 CONTINUE 
RETURN 
END 
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